# -------------------------------------------------------------------
# -------------------------------------------------------------------
# Carga de paquetes necesarios para el trabajo práctico
# Instalar con: install.packages(c("ggplot2","dplyr","signal","tuneR","patchwork"))
# -------------------------------------------------------------------
library(ggplot2)    # Visualización de datos (gramática de gráficos)
library(dplyr)      # Manipulación de data frames
library(signal)     # Procesamiento de señales (FFT, filtros)
install.packages("patchwork")
# -------------------------------------------------------------------
# Configuración global del documento
# Opciones que aplican a todos los chunks de código a continuación
# -------------------------------------------------------------------
# -------------------------------------------------------------------
# Carga de paquetes necesarios para el trabajo práctico
# Instalar con: install.packages(c("ggplot2","dplyr","signal","tuneR","patchwork"))
# -------------------------------------------------------------------
library(ggplot2)    # Visualización de datos (gramática de gráficos)
library(dplyr)      # Manipulación de data frames
library(signal)     # Procesamiento de señales (FFT, filtros)
install.packages("signal")
install.packages("tuneR")
# -------------------------------------------------------------------
# Configuración global del documento
# Opciones que aplican a todos los chunks de código a continuación
# -------------------------------------------------------------------
# -------------------------------------------------------------------
# Carga de paquetes necesarios para el trabajo práctico
# Instalar con: install.packages(c("ggplot2","dplyr","signal","tuneR","patchwork"))
# -------------------------------------------------------------------
library(ggplot2)    # Visualización de datos (gramática de gráficos)
library(dplyr)      # Manipulación de data frames
library(signal)     # Procesamiento de señales (FFT, filtros)
library(patchwork)  # Combinar múltiples gráficos ggplot en una figura
library(scales)     # Formateo de ejes (notación científica, etc.)
# -------------------------------------------------------------------
# FIGURA 2: Transformada de Fourier de las tres señales
# Se usa fft() de R base y se convierte a espectro de potencia en dB.
# -------------------------------------------------------------------
# Función auxiliar: calcula el espectro de potencia unilateral en dB
# Argumento 'x': vector de la señal completa (no recortado)
# Argumento 'fs': frecuencia de muestreo en Hz
# Retorna: data frame con columnas 'freq' (Hz) y 'dB'
calcular_espectro <- function(x, fs) {
N   <- length(x)               # número total de muestras
X   <- fft(x)                  # DFT compleja (N puntos)
# Módulo al cuadrado = potencia; factor 2 por espectro unilateral
Pot <- (2 / N^2) * Mod(X[1:(N/2)])^2
# Convertir a dB (referencia: máximo de la señal)
dB  <- 10 * log10(Pot / max(Pot) + 1e-12)
# Vector de frecuencias correspondiente
freq <- (0:(N/2 - 1)) * fs / N
data.frame(freq = freq, dB = dB)
}
# Calcular espectros usando la señal COMPLETA (no recortada) para mejor resolución
esp1    <- calcular_espectro(s1,     fs) %>% mutate(señal = paste0("s₁ (", f1, " Hz)"))
# Función auxiliar: calcula el espectro de potencia unilateral en dB
# Argumento 'x': vector de la señal completa (no recortado)
# Argumento 'fs': frecuencia de muestreo en Hz
# Retorna: data frame con columnas 'freq' (Hz) y 'dB'
calcular_espectro <- function(x, fs) {
N   <- length(x)               # número total de muestras
X   <- fft(x)                  # DFT compleja (N puntos)
# Módulo al cuadrado = potencia; factor 2 por espectro unilateral
Pot <- (2 / N^2) * Mod(X[1:(N/2)])^2
# Convertir a dB (referencia: máximo de la señal)
dB  <- 10 * log10(Pot / max(Pot) + 1e-12)
# Vector de frecuencias correspondiente
freq <- (0:(N/2 - 1)) * fs / N
data.frame(freq = freq, dB = dB)
}
# Calcular espectros usando la señal COMPLETA (no recortada) para mejor resolución
esp1    <- calcular_espectro(s1,     fs) %>% mutate(señal = paste0("s₁ (", f1, " Hz)"))
# -------------------------------------------------------------------
# FIGURA 2: Transformada de Fourier de las tres señales
# Se usa fft() de R base y se convierte a espectro de potencia en dB.
# -------------------------------------------------------------------
# Función auxiliar: calcula el espectro de potencia unilateral en dB
# Argumento 'x': vector de la señal completa (no recortado)
# Argumento 'fs': frecuencia de muestreo en Hz
# Retorna: data frame con columnas 'freq' (Hz) y 'dB'
calcular_espectro <- function(x, fs) {
N   <- length(x)               # número total de muestras
X   <- fft(x)                  # DFT compleja (N puntos)
# Módulo al cuadrado = potencia; factor 2 por espectro unilateral
Pot <- (2 / N^2) * Mod(X[1:(N/2)])^2
# Convertir a dB (referencia: máximo de la señal)
dB  <- 10 * log10(Pot / max(Pot) + 1e-12)
# Vector de frecuencias correspondiente
freq <- (0:(N/2 - 1)) * fs / N
data.frame(freq = freq, dB = dB)
}
# Calcular espectros usando la señal COMPLETA (no recortada) para mejor resolución
esp1    <- calcular_espectro(s1,     fs) %>% mutate(señal = paste0("s₁ (", f1, " Hz)"))
# -------------------------------------------------------------------
# FIGURA 1: Las tres señales en el dominio temporal
# Se usan tres paneles apilados con patchwork para facilitar la comparación.
# -------------------------------------------------------------------
# Paleta de colores consistente con la convención del marco teórico
col_s1     <- "#1565C0"   # azul — señal 1
col_s2     <- "#C62828"   # rojo — señal 2
col_suma   <- "#2E7D32"   # verde — señal suma
col_env    <- "#7B1FA2"   # violeta — envolvente
# --- Panel superior: Señal 1 ---
p1 <- ggplot(df, aes(x = t)) +
geom_line(aes(y = s1), colour = col_s1, linewidth = 0.6) +
labs(
title = bquote(s[1](t) == A[1] * sin(2*pi * .(f1) * t)),
x = NULL, y = "Amplitud"
) +
scale_x_continuous(labels = scales::number_format(accuracy = 0.001)) +
theme_classic(base_size = 9) +
theme(plot.title = element_text(size = 9, face = "plain"))
# -------------------------------------------------------------------
# FIGURA 1: Las tres señales en el dominio temporal
# Se usan tres paneles apilados con patchwork para facilitar la comparación., fig.cap="Señales individuales $s_1(t)$, $s_2(t)$ y su superposición $s(t) = s_1+s_2$. Se observa cómo la amplitud de la señal suma varía lentamente a la frecuencia de batido.", fig.height=5
# -------------------------------------------------------------------
# Paleta de colores consistente con la convención del marco teórico
col_s1     <- "#1565C0"   # azul — señal 1
col_s2     <- "#C62828"   # rojo — señal 2
col_suma   <- "#2E7D32"   # verde — señal suma
col_env    <- "#7B1FA2"   # violeta — envolvente
# --- Panel superior: Señal 1 ---
p1 <- ggplot(df, aes(x = t)) +
geom_line(aes(y = s1), colour = col_s1, linewidth = 0.6) +
labs(
title = bquote(s[1](t) == A[1] * sin(2*pi * .(f1) * t)),
x = NULL, y = "Amplitud"
) +
scale_x_continuous(labels = scales::number_format(accuracy = 0.001)) +
theme_classic(base_size = 9) +
theme(plot.title = element_text(size = 9, face = "plain"))
library(ggplot2)
# -------------------------------------------------------------------
# FIGURA 1: Las tres señales en el dominio temporal
# Se usan tres paneles apilados con patchwork para facilitar la comparación., fig.cap="Señales individuales $s_1(t)$, $s_2(t)$ y su superposición $s(t) = s_1+s_2$. Se observa cómo la amplitud de la señal suma varía lentamente a la frecuencia de batido.", fig.height=5
# -------------------------------------------------------------------
# Paleta de colores consistente con la convención del marco teórico
col_s1     <- "#1565C0"   # azul — señal 1
col_s2     <- "#C62828"   # rojo — señal 2
col_suma   <- "#2E7D32"   # verde — señal suma
col_env    <- "#7B1FA2"   # violeta — envolvente
# --- Panel superior: Señal 1 ---
p1 <- ggplot(df, aes(x = t)) +
geom_line(aes(y = s1), colour = col_s1, linewidth = 0.6) +
labs(
title = bquote(s[1](t) == A[1] * sin(2*pi * .(f1) * t)),
x = NULL, y = "Amplitud"
) +
scale_x_continuous(labels = scales::number_format(accuracy = 0.001)) +
theme_classic(base_size = 9) +
theme(plot.title = element_text(size = 9, face = "plain"))
df
# -------------------------------------------------------------------
# GENERACIÓN DE LAS SEÑALES EN EL DOMINIO TEMPORAL
# Se crea un vector de tiempo y luego las señales analíticas
# -------------------------------------------------------------------
# Vector de tiempo: desde 0 hasta 'duracion' con paso 1/fs
t <- seq(0, duracion - 1/fs, by = 1/fs)
# -------------------------------------------------------------------
# PARÁMETROS DE LA SIMULACIÓN
# Modificar estos valores para reproducir diferentes escenarios de batido.
# En el laboratorio: reemplazar por los valores medidos con el RTL-SDR.
# -------------------------------------------------------------------
# Frecuencias de las dos señales (en Hz, usando valores de audio para la simulación)
# En el experimento real serían ~144.500 MHz y ~144.501 MHz
f1 <- 440    # Hz — frecuencia de la señal 1 (ej. nota La musical)
f2 <- 443    # Hz — frecuencia de la señal 2 (diferencia = 3 Hz → batido a 3 Hz)
# Amplitudes de cada señal (adimensionales, normalizadas)
A1 <- 1.0   # Amplitud de la señal 1
A2 <- 1.0   # Amplitud de la señal 2
# Parámetros temporales de la simulación
fs   <- 44100   # Frecuencia de muestreo (Hz) — calidad de audio estándar
duracion <- 2   # Duración de la señal simulada (segundos)
# Frecuencias derivadas (calculadas automáticamente)
f_batido     <- abs(f1 - f2)              # Frecuencia de batido perceptible
f_portadora  <- (f1 + f2) / 2            # Frecuencia de la portadora
T_batido     <- 1 / f_batido             # Período del batido (segundos)
# Mostrar resumen de parámetros en la salida del documento
cat(sprintf("f1            = %.1f Hz\n", f1))
cat(sprintf("f2            = %.1f Hz\n", f2))
cat(sprintf("f_batido      = %.1f Hz  (período = %.4f s)\n", f_batido, T_batido))
cat(sprintf("f_portadora   = %.1f Hz\n", f_portadora))
cat(sprintf("Muestras tot. = %d\n", fs * duracion))
# -------------------------------------------------------------------
# GENERACIÓN DE LAS SEÑALES EN EL DOMINIO TEMPORAL
# Se crea un vector de tiempo y luego las señales analíticas
# -------------------------------------------------------------------
# Vector de tiempo: desde 0 hasta 'duracion' con paso 1/fs
t <- seq(0, duracion - 1/fs, by = 1/fs)
# Señal 1: sinusoide pura en f1
s1 <- A1 * sin(2 * pi * f1 * t)
# Señal 2: sinusoide pura en f2
s2 <- A2 * sin(2 * pi * f2 * t)
# Señal suma: superposición lineal de s1 y s2
# Esta es la señal que físicamente llega al receptor SDR
s_suma <- s1 + s2
# Señal de batido analítica (usando la identidad trigonométrica del marco teórico)
# envolvente * portadora — permite graficar la envolvente por separado
s_envolvente <- 2 * A1 * cos(2 * pi * (f_batido / 2) * t)   # envolvente lenta
s_batido_analitico <- s_envolvente * sin(2 * pi * f_portadora * t)
# Organizar en un data frame para facilitar el graficado con ggplot2
# Se limita a los primeros N_mostrar puntos para que los gráficos sean legibles
N_mostrar <- round(fs * 3 * T_batido)   # mostrar 3 períodos de batido completos
N_mostrar <- min(N_mostrar, length(t))  # no exceder el vector real
df <- data.frame(
t          = t[1:N_mostrar],       # eje temporal (s)
s1         = s1[1:N_mostrar],      # señal 1
s2         = s2[1:N_mostrar],      # señal 2
s_suma     = s_suma[1:N_mostrar],  # señal suma
envolvente = s_envolvente[1:N_mostrar]  # envolvente del batido
)
library(ggplot2)
# -------------------------------------------------------------------
# FIGURA 1: Las tres señales en el dominio temporal
# Se usan tres paneles apilados con patchwork para facilitar la comparación., fig.cap="Señales individuales $s_1(t)$, $s_2(t)$ y su superposición $s(t) = s_1+s_2$. Se observa cómo la amplitud de la señal suma varía lentamente a la frecuencia de batido.", fig.height=5
# -------------------------------------------------------------------
# Paleta de colores consistente con la convención del marco teórico
col_s1     <- "#1565C0"   # azul — señal 1
col_s2     <- "#C62828"   # rojo — señal 2
col_suma   <- "#2E7D32"   # verde — señal suma
col_env    <- "#7B1FA2"   # violeta — envolvente
# --- Panel superior: Señal 1 ---
p1 <- ggplot(df, aes(x = t)) +
geom_line(aes(y = s1), colour = col_s1, linewidth = 0.6) +
labs(
title = bquote(s[1](t) == A[1] * sin(2*pi * .(f1) * t)),
x = NULL, y = "Amplitud"
) +
scale_x_continuous(labels = scales::number_format(accuracy = 0.001)) +
theme_classic(base_size = 9) +
theme(plot.title = element_text(size = 9, face = "plain"))
# --- Panel medio: Señal 2 ---
p2 <- ggplot(df, aes(x = t)) +
geom_line(aes(y = s2), colour = col_s2, linewidth = 0.6) +
labs(
title = bquote(s[2](t) == A[2] * sin(2*pi * .(f2) * t)),
x = NULL, y = "Amplitud"
) +
scale_x_continuous(labels = scales::number_format(accuracy = 0.001)) +
theme_classic(base_size = 9) +
theme(plot.title = element_text(size = 9, face = "plain"))
# --- Panel inferior: Señal suma con su envolvente ---
p3 <- ggplot(df, aes(x = t)) +
# Primero la envolvente (positiva y negativa) para que quede detrás
geom_line(aes(y =  envolvente), colour = col_env,
linewidth = 0.8, linetype = "dashed") +
geom_line(aes(y = -envolvente), colour = col_env,
linewidth = 0.8, linetype = "dashed") +
# Luego la señal suma encima
geom_line(aes(y = s_suma), colour = col_suma, linewidth = 0.6) +
labs(
title = bquote(s(t) == s[1] + s[2] ~ "  [envolvente en violeta]"),
x = "Tiempo (s)", y = "Amplitud"
) +
scale_x_continuous(labels = scales::number_format(accuracy = 0.001)) +
theme_classic(base_size = 9) +
theme(plot.title = element_text(size = 9, face = "plain"))
# Combinar los tres paneles verticalmente con patchwork
# plot_layout(ncol=1) apila, heights controla la proporción de cada panel
(p1 / p2 / p3) + plot_layout(ncol = 1, heights = c(1, 1, 1.2))
# -------------------------------------------------------------------
# FIGURA 2: Transformada de Fourier de las tres señales
# Se usa fft() de R base y se convierte a espectro de potencia en dB.
# -------------------------------------------------------------------
# Función auxiliar: calcula el espectro de potencia unilateral en dB
# Argumento 'x': vector de la señal completa (no recortado)
# Argumento 'fs': frecuencia de muestreo en Hz
# Retorna: data frame con columnas 'freq' (Hz) y 'dB'
calcular_espectro <- function(x, fs) {
N   <- length(x)               # número total de muestras
X   <- fft(x)                  # DFT compleja (N puntos)
# Módulo al cuadrado = potencia; factor 2 por espectro unilateral
Pot <- (2 / N^2) * Mod(X[1:(N/2)])^2
# Convertir a dB (referencia: máximo de la señal)
dB  <- 10 * log10(Pot / max(Pot) + 1e-12)
# Vector de frecuencias correspondiente
freq <- (0:(N/2 - 1)) * fs / N
data.frame(freq = freq, dB = dB)
}
# Calcular espectros usando la señal COMPLETA (no recortada) para mejor resolución
esp1    <- calcular_espectro(s1,     fs) %>% mutate(señal = paste0("s₁ (", f1, " Hz)"))
esp2    <- calcular_espectro(s2,     fs) %>% mutate(señal = paste0("s₂ (", f2, " Hz)"))
esp_sum <- calcular_espectro(s_suma, fs) %>% mutate(señal = "Suma s₁+s₂")
# Combinar para graficar con ggplot2 (formato largo)
df_esp <- bind_rows(esp1, esp2, esp_sum)
# Limitar el eje X a una ventana estrecha alrededor de f1 y f2
# para que los picos sean claramente visibles
margen <- max(50, 20 * f_batido)   # margen en Hz
ggplot(df_esp %>% filter(freq >= min(f1,f2) - margen,
freq <= max(f1,f2) + margen),
aes(x = freq, y = dB, colour = señal)) +
# Líneas de referencia verticales en f1 y f2
geom_vline(xintercept = f1, linetype = "dotted", colour = col_s1, alpha = 0.5) +
geom_vline(xintercept = f2, linetype = "dotted", colour = col_s2, alpha = 0.5) +
geom_line(linewidth = 0.7) +
# Etiquetas de frecuencia
annotate("text", x = f1, y = -2, label = paste0("f₁=", f1, "Hz"),
hjust = 1.1, size = 2.5, colour = col_s1) +
annotate("text", x = f2, y = -2, label = paste0("f₂=", f2, "Hz"),
hjust = -0.1, size = 2.5, colour = col_s2) +
scale_colour_manual(
values = c(
paste0("s₁ (", f1, " Hz)") = col_s1,
# -------------------------------------------------------------------
# FIGURA 2: Transformada de Fourier de las tres señales
# Se usa fft() de R base y se convierte a espectro de potencia en dB.
# -------------------------------------------------------------------
# Función auxiliar: calcula el espectro de potencia unilateral en dB
# Argumento 'x': vector de la señal completa (no recortado)
# Argumento 'fs': frecuencia de muestreo en Hz
# Retorna: data frame con columnas 'freq' (Hz) y 'dB'
calcular_espectro <- function(x, fs) {
N   <- length(x)               # número total de muestras
X   <- fft(x)                  # DFT compleja (N puntos)
# Módulo al cuadrado = potencia; factor 2 por espectro unilateral
Pot <- (2 / N^2) * Mod(X[1:(N/2)])^2
# Convertir a dB (referencia: máximo de la señal)
dB  <- 10 * log10(Pot / max(Pot) + 1e-12)
# Vector de frecuencias correspondiente
freq <- (0:(N/2 - 1)) * fs / N
data.frame(freq = freq, dB = dB)
}
# Calcular espectros usando la señal COMPLETA (no recortada) para mejor resolución
esp1    <- calcular_espectro(s1,     fs) %>% mutate(señal = paste0("s₁ (", f1, " Hz)"))
esp2    <- calcular_espectro(s2,     fs) %>% mutate(señal = paste0("s₂ (", f2, " Hz)"))
esp_sum <- calcular_espectro(s_suma, fs) %>% mutate(señal = "Suma s₁+s₂")
# Combinar para graficar con ggplot2 (formato largo)
df_esp <- bind_rows(esp1, esp2, esp_sum)
# Limitar el eje X a una ventana estrecha alrededor de f1 y f2
# para que los picos sean claramente visibles
margen <- max(50, 20 * f_batido)   # margen en Hz
ggplot(df_esp %>% filter(freq >= min(f1,f2) - margen,
freq <= max(f1,f2) + margen),
aes(x = freq, y = dB, colour = señal)) +
# Líneas de referencia verticales en f1 y f2
geom_vline(xintercept = f1, linetype = "dotted", colour = col_s1, alpha = 0.5) +
geom_vline(xintercept = f2, linetype = "dotted", colour = col_s2, alpha = 0.5) +
geom_line(linewidth = 0.7) +
# Etiquetas de frecuencia
annotate("text", x = f1, y = -2, label = paste0("f₁=", f1, "Hz"),
hjust = 1.1, size = 2.5, colour = col_s1) +
annotate("text", x = f2, y = -2, label = paste0("f₂=", f2, "Hz"),
hjust = -0.1, size = 2.5, colour = col_s2) +
scale_colour_manual(
values = c(
paste0("s₁ (", f1, " Hz)") = col_s1,
# -------------------------------------------------------------------
# FIGURA 3: Espectrograma de la señal suma
# Permite visualizar simultáneamente el dominio temporal y frecuencial.
# Se calcula con la STFT (Short-Time Fourier Transform) de forma manual., fig.cap="Espectrograma de la señal suma. El color representa la potencia espectral en cada instante de tiempo. Se observan las dos componentes de frecuencia constante, cuya interferencia genera la modulación de amplitud visible en la Fig.\\ 1.", fig.height=3.5
# -------------------------------------------------------------------
# Parámetros de la STFT
N_ventana <- 2048       # tamaño de la ventana de análisis (muestras)
N_overlap  <- 1024      # solapamiento entre ventanas consecutivas
N_fft      <- 4096      # puntos de FFT (zero-padding para mayor resolución)
hop        <- N_ventana - N_overlap  # paso entre ventanas
# Ventana de Hanning: reduce el efecto de fuga espectral
ventana_hanning <- 0.5 * (1 - cos(2 * pi * (0:(N_ventana-1)) / (N_ventana - 1)))
# Calcular la STFT de la señal suma completa
n_frames <- floor((length(s_suma) - N_ventana) / hop) + 1
freq_vec  <- (0:(N_fft/2 - 1)) * fs / N_fft   # vector de frecuencias
# Pre-alocar la matriz del espectrograma
spec_mat <- matrix(0, nrow = N_fft/2, ncol = n_frames)
for (k in seq_len(n_frames)) {
inicio   <- (k - 1) * hop + 1          # índice de inicio de la ventana
fin      <- inicio + N_ventana - 1      # índice de fin
segmento <- s_suma[inicio:fin] * ventana_hanning  # aplicar ventana
X_k      <- fft(c(segmento, rep(0, N_fft - N_ventana)))  # FFT con zero-padding
spec_mat[, k] <- (2/N_ventana^2) * Mod(X_k[1:(N_fft/2)])^2  # potencia
}
# Convertir a dB
spec_dB <- 10 * log10(spec_mat + 1e-12)
spec_dB_norm <- spec_dB - max(spec_dB)  # normalizar a 0 dB el máximo
# Vector de tiempos central de cada ventana (en segundos)
t_frames <- ((1:n_frames - 1) * hop + N_ventana/2) / fs
# Construir data frame para ggplot2 (formato largo: una fila por píxel)
# Filtrar para mostrar solo la región de interés en frecuencia
margen_esp <- max(100, 30 * f_batido)
mask_freq  <- freq_vec >= (min(f1, f2) - margen_esp) &
freq_vec <= (max(f1, f2) + margen_esp)
# Crear el data frame con expand.grid (producto cartesiano t × f)
df_spec <- expand.grid(
t    = t_frames,
freq = freq_vec[mask_freq]
)
# Rellenar con los valores de potencia correspondientes
spec_subset    <- spec_dB_norm[mask_freq, ]
df_spec$potencia <- as.vector(spec_subset)  # aplanar columna a columna
ggplot(df_spec, aes(x = t, y = freq, fill = potencia)) +
geom_tile() +  # cada "tile" es un píxel del espectrograma
scale_fill_gradientn(
colours = c("navy", "blue", "cyan", "yellow", "red"),
limits  = c(-60, 0),
na.value = "navy",
name    = "dB"
) +
geom_hline(yintercept = c(f1, f2), linetype = "dashed",
colour = "white", alpha = 0.7, linewidth = 0.4) +
annotate("text", x = max(t_frames)*0.02, y = f1,
label = paste0("f₁=", f1, "Hz"), vjust = -0.4,
colour = "white", size = 2.5) +
annotate("text", x = max(t_frames)*0.02, y = f2,
label = paste0("f₂=", f2, "Hz"), vjust = 1.2,
colour = "white", size = 2.5) +
labs(
title = "Espectrograma de la señal suma (STFT)",
x     = "Tiempo (s)",
y     = "Frecuencia (Hz)"
) +
theme_classic(base_size = 9) +
theme(legend.key.height = unit(0.4, "cm"),
legend.text        = element_text(size = 7))
# -------------------------------------------------------------------
# TABLA COMPARATIVA: valores teóricos (calculados con la identidad
# trigonométrica) vs. valores obtenidos por FFT de la señal simulada.
# En el experimento real, la columna "Medido (RTL-SDR)" se rellena
# con los valores leídos en SDR# o en el archivo IQ capturado.
# -------------------------------------------------------------------
# Localizar el pico de mayor potencia en el espectro de la señal suma
# dentro de una ventana de ±50 Hz alrededor de cada frecuencia esperada
encontrar_pico <- function(espectro, f_esperada, margen = 50) {
mask <- espectro$freq >= (f_esperada - margen) &
espectro$freq <= (f_esperada + margen)
espectro$freq[mask][which.max(espectro$dB[mask])]
}
# Espectro de la señal suma completa (recalcular para mayor resolución)
esp_suma_full <- calcular_espectro(s_suma, fs)
f1_medida <- encontrar_pico(esp_suma_full, f1)
f2_medida <- encontrar_pico(esp_suma_full, f2)
tabla <- data.frame(
Parámetro           = c("f₁ (Hz)", "f₂ (Hz)",
"f_batido = |f₁−f₂| (Hz)",
"f_portadora = (f₁+f₂)/2 (Hz)",
"Período batido T = 1/f_batido (ms)"),
`Teórico`            = c(f1, f2,
f_batido,
f_portadora,
round(T_batido * 1000, 2)),
`FFT / Simulado`     = c(f1_medida, f2_medida,
round(abs(f1_medida - f2_medida), 2),
round((f1_medida + f2_medida)/2, 2),
round(1/abs(f1_medida - f2_medida)*1000, 2)),
`Medido (RTL-SDR)`   = c("___", "___", "___", "___", "___"),
check.names = FALSE
)
knitr::kable(tabla,
caption = "Comparación de parámetros teóricos, simulados y experimentales.",
align   = c("l", "r", "r", "r"),
booktabs = TRUE   # mejor estilo para PDF LaTeX
)
# -------------------------------------------------------------------
# ERROR RELATIVO entre el valor teórico y el obtenido por FFT.
# Completar la columna "Medido" con el valor real del RTL-SDR.
# -------------------------------------------------------------------
# Ejemplo de cálculo del error para la frecuencia de batido
f_bat_fft <- abs(f1_medida - f2_medida)
error_pct <- abs(f_bat_fft - f_batido) / f_batido * 100
cat(sprintf("Frecuencia de batido teórica  : %.4f Hz\n", f_batido))
cat(sprintf("Frecuencia de batido por FFT  : %.4f Hz\n", f_bat_fft))
cat(sprintf("Error relativo                : %.4f %%\n", error_pct))
# El error es función de la resolución espectral Δf = fs / N
resolucion_FFT <- fs / length(s_suma)
cat(sprintf("Resolución espectral (Δf)     : %.4f Hz\n", resolucion_FFT))
# -------------------------------------------------------------------
# Información de la sesión de R (reproducibilidad)
# -------------------------------------------------------------------
cat("Versión de R:", R.version.string, "\n")
cat("Paquetes clave:\n")
cat(sprintf("  ggplot2  %s\n", packageVersion("ggplot2")))
cat(sprintf("  signal   %s\n", packageVersion("signal")))
cat(sprintf("  patchwork %s\n", packageVersion("patchwork")))
# Función auxiliar: calcula el espectro de potencia unilateral en dB
# Argumento 'x': vector de la señal completa (no recortado)
# Argumento 'fs': frecuencia de muestreo en Hz
# Retorna: data frame con columnas 'freq' (Hz) y 'dB'
calcular_espectro <- function(x, fs) {
N   <- length(x)               # número total de muestras
X   <- fft(x)                  # DFT compleja (N puntos)
# Módulo al cuadrado = potencia; factor 2 por espectro unilateral
Pot <- (2 / N^2) * Mod(X[1:(N/2)])^2
# Convertir a dB (referencia: máximo de la señal)
dB  <- 10 * log10(Pot / max(Pot) + 1e-12)
# Vector de frecuencias correspondiente
freq <- (0:(N/2 - 1)) * fs / N
data.frame(freq = freq, dB = dB)
}
# Calcular espectros usando la señal COMPLETA (no recortada) para mejor resolución
esp1    <- calcular_espectro(s1,     fs) %>% mutate(señal = paste0("s₁ (", f1, " Hz)"))
esp2    <- calcular_espectro(s2,     fs) %>% mutate(señal = paste0("s₂ (", f2, " Hz)"))
esp_sum <- calcular_espectro(s_suma, fs) %>% mutate(señal = "Suma s₁+s₂")
# Combinar para graficar con ggplot2 (formato largo)
df_esp <- bind_rows(esp1, esp2, esp_sum)
# Limitar el eje X a una ventana estrecha alrededor de f1 y f2
# para que los picos sean claramente visibles
margen <- max(50, 20 * f_batido)   # margen en Hz
ggplot(df_esp %>% filter(freq >= min(f1,f2) - margen,
freq <= max(f1,f2) + margen),
aes(x = freq, y = dB, colour = señal)) +
# Líneas de referencia verticales en f1 y f2
geom_vline(xintercept = f1, linetype = "dotted", colour = col_s1, alpha = 0.5) +
geom_vline(xintercept = f2, linetype = "dotted", colour = col_s2, alpha = 0.5) +
geom_line(linewidth = 0.7) +
# Etiquetas de frecuencia
annotate("text", x = f1, y = -2, label = paste0("f₁=", f1, "Hz"),
hjust = 1.1, size = 2.5, colour = col_s1) +
annotate("text", x = f2, y = -2, label = paste0("f₂=", f2, "Hz"),
hjust = -0.1, size = 2.5, colour = col_s2) +
scale_colour_manual(
values = c(
paste0("s₁ (", f1, " Hz)") = col_s1,
df_esp
